Quantum Monte Carlo simulations of bosonic and fermionic impurities 
in a two-dimensional hard-core boson system 
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A two-dimensional lattice hard-core boson system with a small fraction of bosonic or fermionic 
impurity particles is studied. The impurities have the same hopping and interactions as the dominant 
bosons and their effects are solely due to quantum statistics. Quantum Monte Carlo simulations are 
carried out in which paths of the dominant boson species are sampled and a summation is performed 
over all second-species paths compatible with the permutation cycles. Both kinds of impurities 
reduce modestly and equally the Kosterliz-Thouless superfluid transition temperature. However, 
the effective impurity interactions are found to be qualitatively different at lower temperatures; 
fermions are repulsive and further suppress superfluidity as T — > 0. 
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Ultracold atoms in optical lattices offer unprecedented 
opportunities to realize novel quantum states of matter 
[l[ . One prospect is to tailor systems to mimic hamiltoni- 
ans of fundamental interest in condensed matter physics, 
e.g., the Hubbard model [2j|. Another interesting route is 
to explore states that do not have any realizations in nat- 
urally occurring systems. The use of mixtures of different 
atomic species open up almost endless possibilities. In 
the case of two species, there are three mixture classes; 
bose-bose, fermi-fermi, and bose-fermi, with the latter 
perhaps offering the most interesting prospects. Several 
exotic phases have been predicted theoretically, e.g., su- 
persolids [§| , several different Mott states [1] , multiply de- 
generate quantum-disordered states with glass-like prop- 
erties [Bj], paired states with various orbital symmetry 
Q , and a host of states of effective fermion-boson com- 
posite particles 0]. Experimentally, there are intriguing 
results indicating a strong influence of a small admixture 
of fermions in boson gases in optical lattices [H . 

In this Letter a simple two-species model will be con- 
sidered in which all particles have identical hoppings and 
interactions, posing a clean way to elucidate the funda- 
mental role of quantum statistics. The situation can be 
realized experimentally in systems with two isotopes of 
the same atom, e.g., 6 Li- 7 Li mixtures A ID model 
of this kind has been studied using various analytical 
approaches [lj|. Here a quantum Monte Carlo (QMC) 
method is developed for a low concentration of fermionic 
or bosonic impurities in a bath of bosons. The scheme 
is applied to a 2D model. Defining creation operators a\ 
and b\ for the two species the hamiltonian is 



<ii> 



+ a]a l + b]b j + b]b i ), 



(1) 



where (ij) denotes nearest neighbors on a square lattice 
of N = L x L sites with periodic boundaries. Species A is 
bosonic, whereas the B particles can be either bosons or 
fermions. Both are subject to a hard-core constraint, i.e., 
the sites are either empty or singly-occupied. Half-filling 



in the canonical ensemble will be considered here; n = 
riA + n B = N/2, where ua = a l a i an d n B — 
The trapping potential necessary to model experiments 
is left out, but can easily be incorporated in future work 
using the QMC scheme introduced below. 

With a fermionic B species, the hamiltonian as a func- 
tion of ris /n interpolates between spinless fermions and 
the standard hard-core boson model, which is equivalent 
to the S — 1/2 XY model and undergoes a Kosterliz- 
Thouless (KT) transition to a superfluid with power-law 
off-diagonal correlations at temperature TktA ~ 0.69 
ll|: at T = it is long-range ordered. The spinless 



fermion model, on the other hand, has a metallic ground 
state and does not undergo any finite-T transition. It is 
then interesting to consider the evolution of the ground 
state and finite-T properties as a function of the fermion 
fraction. Here the A species will be considered dominant, 
the B particles acting as impurities; jjb < ua- 

The purely bosonic ground state does not change when 
some A particles are replaced by another boson species. 
The excitations are affected, however, and one can ex- 
pect a reduction in Tkt relative to the single-species 
system. It will be shown here that the effective inter- 
action between the impurities is attractive at high tem- 
peratures but changes, in a singular way, to repulsive at 
Tkt- Bosonic impurities become attractive at lower T, 
whereas fermions stay repulsive as T — > 0. The effects of 
the impurities on the KT transition are independent of 
their statistics, but there are indications of a more dra- 
matic suppression of superfluidity by fermions as T — > 0. 

QMC Algorithm— Consider finite-T QMC methods in 
which the density matrix e~ H / T (kg = 1) is written as 
a sum of operators Pj which propagate real-space states 
| a) such that Pi\a) — Wi(a)\a). The paths (i,a) are 
importance-sampled according to their weights Wi{a) 
in the partition function Z = a Wi{a) . Such path- 
integral methods based on "time-slicing" can be formu- 
lated in the continuum [l2| and on lattices [l3| • In recent 
years very efficient algorithms for updating the paths 
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FIG. 1: (Color online) Permutation cycles for an 8 x 8 sys- 
tem at T/t = 0.5. Open and solid circles represent empty 
sites and particles, respectively. The net effect of the SSE 
propagation of the state is to cyclically permute the particles 
connected by the closed paths. The isolated solid (black) cir- 
cles are particles not undergoing any net permutations; they 
are considered permutation cycles of length one. 



have been devised [1J, [15|, [16J and systems with thou- 



sands of bosons can now be routinely simulated even at 
low temperatures. On the lattice, the alternative stochas- 
tic series expansion (SSE) approach [IE> lU 13 > in which 
Q-H/T j g Baylor expanded to all contributing orders, is of- 
ten more efficient and will be employed here. The details 
of the method are unimportant, however, and the scheme 
for treating impurity particles outlined below should be 
applicable with any standard boson QMC algorithm. 

Identical particles in Pi\a) are permutations of those in 
| a). The permutations can be decomposed into cycles Cj 
in which rrii particles are permuted independently of the 
other particles. An example from an actual SSE simula- 
tion of an 8 x 8 lattice with 32 bosons is shown in Fig. [1] 
In the case of a multi-species system, a term (path) con- 
tributes to the partition function as long as all particles 
within each individual permutation cycle are identical; 
particles in different cycles do not have to be identical. 
Now, if all the interaction and hopping parameters are 
identical for all species, the weight Wi(a) is independent 
of the (allowed) distributions of the particles over the 
cycles. A simulation can then be carried out for all iden- 
tical hard-core bosons (A particles) and, subsequently, 
when measuring observables, two (or more) species can 
be considered by substituting all the particles in some of 
the cycles with B particles. All ways of filling the cycles 
can be summed up exactly, including fermionic signs. 

Denote by N c (m) the number of cycles of length m and 
by n c (m) the number of cycles filled with B particles. 
Then the total number of ways of distributing rib bosons 
of type B in the system (substituting the same number 
of A particles) is given by 



w b 



En 

{n c } m=l 



N c (m) 
n c (m) 



(2) 



where, for a fixed number rib of B particles, the 
sum is over all cycle fillings satisfying the constraint 
^2 m mn c (m) = rib. Thus, if the paths are importance- 




FIG. 2: (Color online) (a) Weight distribution for fermions 
(wf > and uif < separately) and bosons in simulations of 
an L = 16 lattice with nij = 12 at T/t = 0.25. Bin k (inte- 
ger) represents the probability of k < log 2 (|ui|) < k + 1 (the 
special case w = is in the k = — 1 bin), (b) The probability 
times the weight, i.e., the actual contribution of terms in a 
given weight range (normalized to unity for both bosons and 
fermions). The scheme breaks down when a significant frac- 
tion of the weight in a histogram in (b) extends far into the 
rarely sampled right tail of the probability histogram in (a). 
Here the average fermionic sign, i.e., the ratio of fermionic 
and bosonic weights, (S) — {wf)/{wb} ~ 0.03. The effective 
sign in the simulation is the difference between the Wf > 
and Wf < histograms in (b) and is larger; (S%ff ~ 0.4 

sampled using the weight Wi(a), the measurements of 
some observable O should be weighted by Wb] 



(0) = 



(3) 



Here the sums are over the paths actually sampled in the 
simulation. The estimator {Oi yCt ) is an average over the 
w b (i,a) different cycle fillings. 

In the case of fermions, each cyclic permutation of an 
even number of particles yields a minus sign. Thus, in 
the presence of n/ fermionic B particles the weight is 



Wf 



En 

{n c } m=l 



N c (m) 
n c (m) 



(-1) 



(m— l)n c (m) 



(4) 



with the constraint ^ m mn c (m) = n.j. 

It will be demonstrated that it is feasible to evaluate 
exactly the weights ©,(21) even for a relative large num- 
ber of impurity particles — results will be presented for 
n bj < 32 on a 32 x 32 lattice. For some observables the 
estimator Oi tCC is independent of the cycle filling whereas 
in other cases the evaluation of its average may be more 
complicated. In some cases, it may be necessarily to 
carry out a separate sampling of the average estimator 
for each path. Here only the internal energy and the 
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FIG. 3: (Color online) QMC and Lanczos results for the in- 
ternal energy of a 4 x 4 system with n/ = 1,2,3 fermions, 
relative to the energy of the purely bosonic system. 

phase stiffness will be considered, both of which have 
SSE estimators involving operator counts in Pj, which 
are independent of the distribution of A and B particles. 
The energy E = (H) is given simply by the average order 
p of the SSE Taylor expansion; E = -(p)/f3 [17]. The 
stiffness is the second derivative of the energy E(<p) with 
respect to a twist <f> in the boundary condition. It is ob- 
tained by averaging the squared winding number [l^.[l9|. 

A limitation of the cycle summation approach is that 
the relative fluctuations in the weighs Wbj grow as ribj 
increases. Then only a decreasing subset of the gener- 
ated configurations will contribute significantly to com- 
puted quantities. In addition, for fermions there is still 
a sign problem, although the summation over many cy- 
cle fillings with different signs does alleviate it signifi- 
cantly. These issues are illustrated and further discussed 
in Fig. [2] Another problem is that it becomes pro- 
hibitively time consuming to exactly evaluate Wbj for 
large ribj. For moderate ribj, it is possible to obtain 
results on large lattices in temperatures regimes where 
interesting physics takes place, as will be shown next. 

Results. — The correctness of the QMC scheme was con- 
firmed by exact diagonalization results for a 2 x 4 system 
at finite T. Even using momentum conservation, it is dif- 
ficult to completely diagonalize a 4 x 4 system with two 
particle species. The ground state can be obtained using 
the Lanczos method, however. 

An interesting quantity is the change A E{n b j) in 
the internal energy relative to the nbj — energy; 
A E = E(n bJ ) - E(0). Fig.[3]shows how QMC results for 
AE(jif)/rif approach corresponding Lanczos results for 
Uf = 1,2,3. Note that with n/ = 1, fermion anticom- 
mutation does not come into play and the same bosonic 
ground state as with n/ = is obtained as T — > 0. 

An effective interaction energy between impurities 
can be defined by subtracting the single-impurity result 
A E (1); E- mt = A E (n b j)/n b j - A E (1). This quantity is 
shown in Fig. [4] for 8 bosonic and fcrmionic impurities in 
lattices of size L = 8, 16, and 32. At high temperatures 
Ejnt < 0, indicating effectively attractive interactions for 
both bosons and fermions. A singular behavior appears 
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FIG. 4: (Color online) Interaction energy for 8 impurities in 
L — 8, 16, and 32 lattices. The dashed lines indicate Tkt- 



to develop at T ss Tkt, where for large systems both the 
fcrmionic and bosonic interactions first become increas- 
ingly attractive and then suddenly turn repulsive. At 
lower T the fermionic interactions stay repulsive while 
the bosonic ones become attractive again. The behavior 
is very similar for all impurity numbers nf > 2 studied. 
In the case of fermions (bosons) E- mt at low T increases as 
a function of n/ (decreases as a function of rib), as would 
be expected for repulsive (attractive) interactions. When 
negative, Ei nt is seen to decrease as a function of the sys- 
tem size, indicating that the impurities do not demix but 
remain distributed throughout the system. 

Note that n/ is fixed in Fig. |4j i.e., the impurity con- 
centration decreases with increasing L. If the impurities 
do not segregate, it is clear that anticommutation will not 
be significant down to some concentration dependent T 
at which fermions, at their typical separation, can begin 
to permute. Thus, the behavior around Tkt for a low 
impurity density should be a bosonic feature. The qual- 
itative difference between bosonic and fcrmionic impu- 
rities emerges at lower temperature, where the fermions 
remain repulsive down to T — > whereas the bosons be- 
come attractive [and eventually the effective interactions 
vanish as E(nb > 0) — > E(0) for bosons]. 

It appears likely that the singularity seen developing in 
Fig.|4]should move exactly to Tkt as L — > oo and, hence, 
that vortices play a decisive role in the effective interac- 
tions. One can speculate that the each impurity particle 
associates with a vortex. The low-T repulsive fermionic 
interactions may then point to vortex-antivortex pairs 
that increase in size and do not annihilate. There would 
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then be n//2 vortex-antivortex pairs left at T = 0. 

The influence of the bosonic and fermionic impurities 
on the phase stiffness (the superfluid density [3]) around 
Tkt is almost indistinguishable for large system sizes, as 
shown for L = 32 in Fig. 0Ja). In both cases the stiff- 
ness is mildly suppressed, pointing to a modest reduction 
in Tkt for impurity concentrations 1/32 (16 impurities) 
and 1/16 (32 impurities). However, as shown in Fig.^b), 
there are significant differences in the case of two impu- 
rities in a 4 x 4 lattice at low temperatures. As expected, 
for bosonic impurities the stiffness approaches that of 
the single-species system. Two fermions, on the other 
hand, strongly suppress the stiffness as T — > 0. In fact, 
p s — > — co as T — > 0, which can be traced to the mini- 
mum in the energy E(<f>) as a function of the boundary 
twist 4> being slightly away from = 0. In addition to 
this anomaly, the ground state does not have momentum 
q = but is degenerate with q = (±7r/2, 0), (0, ±7r/2). 
These features suggest that the ground state is frus- 
trated, which may again be indicative of the two fermions 
being tied to a vortex-antivortex pair. A drop in the stiff- 
ness is also seen for a small (2 — 4) number of fermions in 
larger lattices, at a T which decreases as L is increased 
for fixed n/. The ground state may thus be insulating. 

Summary and Conclusions. — The QMC method devel- 
oped here enables studies of bosonic as well as fermionic 
impurity particles in a bath of hard-core bosons. Here 
a system with no interactions apart from the hard-core 
constraint was studied, but the method is applicable also 
in the presence of other interactions as long as they are 
equal for all species (inter- and intra-species) . 

In the model studied here, signs of a singular effective 
impurity interaction are seen at Tkt, pointing to a de- 
cisive role of vortices. At low T the effective fermionic 
interactions are repulsive whereas the bosonic ones are 
attractive. A scenario suggested by these findings is that 
impurities associate with vortices. An effectively repul- 
sive fermionic impurity interaction, along with a drop 
observed in the phase stiffness, would then imply that 
some vortex-antivortex pairs remain as T — > 0, poten- 
tially leading to an insulating ground state. For bosonic 
impurities, the ground state is the same as without impu- 
rities. Future studies will address correlations in systems 
with a very small number of impurities at lower temper- 
atures. This should give further insights into the nature 
of the ground state of the bose-fermi mixture. 

In the presence of a trapping potential, the repulsive 
fermions should be expelled to the condensate boundary 
at low temperatures. However, simulations of the grand- 
canonical ensemble, but with fixed n/, demonstrate an 
effectively attractive fermi-bose interactions (the filling 
(n) becomes larger than 1/2). Thus the boundary should 
be in an interesting non-trivial mixed state. 

I would like to thank Anatoli Polkovnikov and Asle 
Sudb0 for useful discussions. This work was supported 
by the NSF under grant No. DMR-0513930. 
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FIG. 5: (Color online) Phase stiffness for (a) L = 4 systems 
with nbj = 0, 2 and (b) L — 32 systems with nbj — 0, 16, 32. 
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